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We present large scale computer simulations of the nonlinear bulk rheology of lamellar phases (smectic liquid crystals) at moder- 
ate to large values of the shear rate (Peclet numbers 10 — 100), in both two and three dimensions. In two dimensions we find that 
modest shear rates align the system and stabilise an almost regular lamellar phase, but high shear rates induce the nucleation and 
proliferation of defects, which in steady state is balanced by the annihilation of defects of opposite sign. The shear rate % at onset 
^vq of this second regime is controlled by thermodynamic and kinetic parameters; we offer a scaling analysis that relates % to a crit- 
ical "capillary number" involving those variables. Within the defect proliferation regime, the defects may be partially annealed 
d by slowly decreasing the applied shear rate; this causes marked memory effects, and history-dependent rheology. Simulations in 
three dimensions show instead shear-induced ordering even at the highest shear rates studied here. This suggests that % shifts 
markedly upward on increasing dimensionality. This may in part reflect the reduced constraints on defect motion, allowing them 
to find and annihilate each other more easily. Residual edge defects in the 3D aligned state mostly point along the flow velocity, 
an orientation impossible in two dimensions. 
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1 Introduction 

Instances of soft complex fluids abound in our everyday lives: 
examples include paint (a colloidal suspension), mayonnaise 
(an emulsion), and engine oil (a polymer solution). Many soft 
materials, including the amphiphiles used in skin creams, fab- 
ric softeners, and some pharmaceutical products, form lamel- 
lar mesophases (also known as smectic liquid crystals) whose 
perfectly ordered state consists of a regular stack of flat lay- 
ers. As well as these lyotropic liquid crystals arising in am- 
phiphilic solutions, thermotropic smectics commonly arise in 
systems of rodlike, mesogenic molecules even without sol- 
vent ' . Such phases are of wide technological importance in 
devices such as electronic displays and sensors. 

The flow behaviour (rheology) of smectic liquid crystals 
affects many industrial and manufacturing processes; con- 
versely, rheological measurements provide important probes 
of their macroscopic material properties. For a simple fluid 
the response to steady applied stress is quantified in the linear 
(Newtonian) regime by a single viscosity. Lamellar phases 
are anisotropic, and thus in principle have several viscosities 
as well as an elastic modulus related to layer compression. 
However, these linear response parameters become somewhat 
irrelevant when defects are present: unless their density is very 
low, the defects soon dominate the stress. Moreover, since the 
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defect pattern is strongly dependent on flow rate, the rheology 
is generically nonlinear. On shearing, the effective viscosity 
(ratio of shear stress <3xy to strain rate y) may either increase 
or decrease, corresponding to shear thickening or shear thin- 
ning behavioui^^l In addition, if a defect texture is present in 
the system at rest, this may show a viscosity divergence at 
small Y — i.e., a yield stress. Finally, in common with other 
complex fluids, lamellar phases can undergo transitions into 
heterogeneous states, for instance to form shear bands (layers 
of unequal (5xy or y)3. 

Lamellar Uquid crystals provide an important paradigm for 
non-Newtonian fluids in which the microstructure of the ma- 
terial is strongly coupled to the flow. Crucially, the defect 
texture may depend on past history - for instance, the strain 
rate in any prior flow that was suddenly switched off. This 
rheological hysteresis means that the precise flow protocol, 
including any pre-shearing applied prior to starting the ex- 
periment, needs to be specified. Such hysteresis can arise in 
polydomain liquid crystals of various types (not just smec- 
tics) such as nematics in which the existence and evolution 
of regions with misaligned directors creates rheological non- 
linearity. A closely related example is the so-called 'region I 
flow' in polymeric nematic liquid crystals, where it is thought 
that a defect network arises, with enhanced viscosity, but with 
marked shear thinning as this structure is sheared out at high 
flow rates ^. 

An intriguing complication in addressing the dynamics of 
layered systems is permeation, where individual molecules 
in the fluid translate while keeping the layerwise order in- 
tact. In steady state, permeation leads to increased dissipa- 
tion and hence to a large viscosity in any geometry where 
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flow necessitates such motion' '. (A well-known example of 
permeation flow is found in cholesterics which are sheared or 
pushed along the direction of the cholesteric helix.) Perme- 
ation is often involved in defect motion such as the growth 
or shrinkage of incomplete layers at the sites of edge disloca- 
tions, and may also be implicated in other forms of topolog- 
ical reconstruction. Such reconstruction arises under shear, 
sometimes to an extreme degree, for instance in the shear- 
induced transition from an ordered lamellar stack to an oruon 
phase in lyotropic smectics . (The onion phase comprises a 
space-filling packing of multilayer vesicles; under some con- 
ditions similar defect textures can arise without 
The mechanism underlying such global reconstructions under 
shear of the layer organization remains unclear Most theo- 
retical predictions have focussed on buckling and other in- 
stabilities of perfectly ordered smectics as the initial trigger 
for onion formation ^^ '^ but the pat hway may also strongly 
involve the physics of defect motion'EIIS. Although we do 
not see onions in the simulations reported here, in two dimen- 
sions we do see a transition from a well ordered phase to a 
defect-rich phase, whose possible link to onion formation de- 
serves to be explored in future work. Another interesting av- 
enue connects with the observation, in block copolymers and 
some lyotropic surfactant systems, of a shear-induced phase 
transitions from a homogeneous state into ordered phases like 
the lamellar oneE^E] However, theory suggests that thermal 
noise (causing local fluctuations about the mean state of order) 
plays a central role in this transition*^. In the present work we 
focus on the simplest case in which such thermal fluctuations 
are ignored, and the time evolution of the microstructure is de- 
terministic (though certainly chaotic and unpredictable once 
defects are present). 

Despite a wealth of experimental results, and the impor- 
tance of flow behaviour to many commercial applications, 
mechanistic understanding of the rheology of lamellar phases 
remains limited. This is largely because of the nonlinear na- 
ture of the equations of motion, and the coupling of the flow to 
different microstructures. These complications either require 
drastic approximations or, as pursued here, use of computer 
simulation. However, once defects are present, such simu- 
lations are computationally expensive. While sheared smec- 
tics have been studied extensively via simulations in the lit- 
eratureESESI computational constraints have meant almost all 
these results are limited to two dimensions. Even in 2D, the 
number of layers simulated has often not been large enough to 
meaningfully explore the dynamics of defect textures which 
can have a controlling influence on rheology. Moreover, for 
technical reasons a number of these studies have used bound- 
ary conditions (BCs) that are not optimal if one wants to pre- 
dict the behaviour of bulk samples, for example using solid 
moving walls to drive the fluid flow rather than periodic BCs. 
This choice does not minimize finite size effects and, when 



combined with lattice Boltzmann algorithms, severely limits 
the shear rates attainable in large systems. Moreover, solid 
walls for the fluid velocity have sometimes been used in con- 
junction with periodic BCs for the order parameter - a 
somewhat inconsistent combination that makes it difficult to 
disentangle bulk and boundary effects. 

In this work, we try to overcome these technical obstacles, 
presenting systematic results for large two dimensional simu- 
lations with fully periodic boundary conditions. Guided by 
those results, we additionally report some very large three 
dimensional simulations at selected points in the parameter 
space. These 3D simulations are, we believe, of unprece- 
dented scale for lamellar phases, although even larger ones 
have recen tly be en used to address the rheology of cubic liq- 
uid crystalsl^SlSSI, 

Below, we describe smectics via a widely-used free energy 
functional '^ ^ 2 3 . 3 .1 1 }\ xhis provides a robust generic framework 
for addressing lamellar ordering without being tied to any spe- 
cific class of materials (although we do become more specific 
when estimating suitable parameters for the simulations). The 
free energy is formulated in terms of a scalar order parameter 
(|)(r). This can be viewed as describing a composition variable 
in a binary fluid system stabilized by added amphiphiles; one 
can then speak of alternating layers of fluid, with interfaces 
between these layers. A similar interpretation applies to block 
copolymers where the composition variable describes the rel- 
ative local density of the two blocks. However, the same order 
parameter could equally describe the local number density of 
the midpoints of a rodlike mesogenic molecule (measured rel- 
ative to the mean) which can develop a periodic density wave 
that represents a thermotropic smectic state. 

The governing equations in our model are a convection- 
diffusion equation for the scalar order parameter, and the 
Navier-Stokes equation describing the hydrodynamics of fluid 
flow. The stress tensor for the fluid contains not only a New- 
tonian viscous term but extra contributions arising from local 
deviations of the order parameter from its state of least free en- 
ergy. Diffusive relaxation of the order parameter then causes 
additional dissipation, resulting in an effective viscosity that 
can be much higher than the intrinsic Newtonian contribu- 
tion. Our dynamical equations are solved using a combina- 
tion of a finite difference scheme for the convection-diffusion 
equation, and the lattice Boltzmann method (LBM)^^ for 
the Navier-Stokes equation^- We use sliding periodic 

(Lees-Edwards) boundary conditions for both the convection- 
diffusion and the Navier-Stokes equations. Note that while 
these boundary conditions impose a fixed integrated strain rate 
across the whole system they do not impose a uniform strain 
rate locally; shear bands can therefore arise if the system so 
chooses. (As it happens, we do not observe shear bands in the 
range of shear rates we address here.) 

Our main finding is that, in two dimensions at least, sheared 
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smectics can show two completely different types of Theologi- 
cal response. At modest rates of shear, where thermodynamic 
forces are relatively strong compared with flow effects, shear 
quickly aligns the system into the perfect lamellar phase pre- 
ferred by thermodynamics, with only a few defects. The rhe- 
ology in this regime is quasi-Newtonian, with shear stress lin- 
ear in the applied shear rate. (There is however a large first 
normal stress difference.) At higher shear rates, larger ther- 
modynamic forces are generated as layers fail to keep up with 
the imposed flow. This leads to the nucleation of many new 
defects, balanced in steady state by annihilation when defects 
of opposite sign encounter each other. 

The defects do not appear as a direct result of spontaneous 
thermal fluctuations. Instead they are a result of an undulation 
instability in our (deterministic) nonlinear dynamics. This oc- 
curs when viscous forces begin to dominate interfacial forces 
causing layers to pinch off and break, creating new defects 
in the process. The onset of a nontrivial topological pattern 
under deterministic evolution equations in a complex fluid is 
similar to the case of sheared binary fluids undergoing phase 
separation . 

The strain rate needed to enter this nonlinear flow regime 
depends on thermodynamic and kinetic parameters. A com- 
prehensive exploration of parameter space is not possible, but 
we give a scaling argument which relates % to those parame- 
ters through an effective capillary number, formed as the ratio 
of viscous to thermodynamic stress. Within the defect prolif- 
eration regime, the macroscopic rheology is generically non- 
hnear, and shows shear thinning. The shear stress appears to 
vary as a sublinear power law, Oxy °= 7"', with m ~ 0.87 for 
typical parameter values, although the dependence of m on 
the parameters precludes interpretation of this as a universal 
scaling law. Such a behaviour is qualitatively reminiscent of 
the one found experimentally in^^, where the numerical value 
of m was however different. In this strongly non-Newtonian 
regime we also find that the rheology is history dependent. 

Finally, and intriguingly, we find that dimensionality ap- 
pears to play a key role. Within our 3D simulations, even 
when performed with parameter values identical to those used 
in 2D, we have so far not been able to enter the nonlinear 
flow regime associated with defect proliferation. Instead, we 
always observe shear-induced ordering into a lamellar phase 
with quasi-Newtonian rheology. Further research, requiring 
careful and systematic parameter steering, will be needed to 
map out quantitatively the onset of the defect proliferation 
regime, if one exists, in three dimensions. 

This paper is organized as follows. In Section |2] we intro- 
duce the model, including its equilibrium properties, the full 
stress tensor and the dynamical equations. We briefly sketch 
how we derive structural information from the order parame- 
ter and describe a method for quantitative analysis of the de- 
fect density. Section [3] describes the mapping between our 



simulation parameters and those for physical systems of inter- 
est. We present results for the two-dimensional case in Section 
4. 1 discussing the morphology as a function of flow rate. We 
also consider average and local shear stress densities and first 
normal stress differences, and explore their relation to the de- 
fect structure. In addition we report Fourier-space information 



(structure factors) under flow. In Section 4.2 we present and 
compare our results for 3D simulations which, by computa- 
tional necessity, are less comprehensive than in 2D. A sum- 
mary with conclusions closes the paper in Section|5] 

2 Model and Definitions 

The equilibrium phase in a system of volume V is found for 
an order parameter field (|)(r) = that minimises the Landau- 
Brazovskii free energy functionalEEHI 



(1) 



In minimizing this free ene rgy at equilibrium, we ignore ther- 
mal fluctuations (treated inE2l32l)^ thereby ad opting a mean- 
field approximation as in previous literaturel^^H^HEII. This 
should be a good approximation far from the ordering tran- 
sition at which (|) first becomes nonuniform. The uniform state 
has (|) = and hence !F =0. 

The parameter range relevant to lamellar phases has K < 
and c > 0. Then, if we suppress the quartic term, the free en- 
ergy density of a sinusoidal lamellar state (|)(r) = (|)|(.cos(k.r) 
obeys f = {a + Kk^ + c^^)(|)|/4. This is minimized by choos- 
ing k = q = ■\/|k|/2c; the free energy density is then / = (a — 
K^/4c)(|)^/4 so that ordering will occur for a < Qc = K^/4c. 
The amplitude is then unbounded without the (positive) 
quartic term, since / can in that manner become infinitely 
negative. Restoring the quartic term now causes saturation 
at (|)^ ^ •\/4(flc — a) /3b, creating also nonlinear harmonics 
in the density profile. The final free energy density obeys 
/ ~ {cLc — af-jb. By similar arguments, the characteristic elas- 
tic modulus associated with microstructural deformations is 
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which is the energy density of the gradient terms only. The 
terms involving ^ and (|)^ in Eq.[T]are unaltered by affine de- 
formations which conserve density locally; hence those terms 
do not contribute to g. 

The dynamics of the order parameter is governed by an 
advection-diffusion equation with a diffusive current of (|) de- 
pendent on the chemical potential /j, the functional derivative 
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of J with respect to (|), as follows: 

+ V(m(^) MV^ ^ = MV^^. (3) 

0(|) 

Here M is a mobility coefficient. A wavevector-dependent 
hnearized collective diffusivity D]. can then be defined via 
isfk = ^k^Dj^i^i;, whose value at the ordering wavevector obeys 
— M{a — ac)/2 (this is negative, signifying the instabil- 
ity). In contrast, an estimate of the molecular diffusivity can 
be made by considering zero wavevector, setting b = —a (see 
below) and examining the relaxation equation for small de- 
viations around the uniformly ordered state at (|) = ±1. This 
gives (j) = 2bMV^(^ so that the molecular diffusivity is of order 
D = 2Mb. 

Treating the fluid as incompressible, the evolution of the 
velocity u obeys the Navier-Stokes equation 

p{d,Ua + (MpVp)Mc() = DpOdp, (4) 

where p is the mass density. The stress tensor in this equation 
reads 

Oap =1l(3aMp+3pMa)-/:'5ap-Pap- (5) 

This consists of a Newtonian hydrodynamic part, which is pro- 
portional to the background viscosity r| and the symmetrized 
velocity gradient, a hydrostatic pressure term, which includes 
an isotropic stress contribution from the order parameter 

+c (^(^(V2)2^ + ^ (V»2 + d^dyiV^)^ (6) 

and finally another term consisting of the o rder-parameter con- 
tribution to the deviatoric pressure^ISlSI 

^ap = KaaH^-c(3aH(^^*)+^P^^«(V>)) . (7) 

In our hybrid LB scheme, the Navier-Stokes equation (j4]) is 
solved at every timestep by means of a multi-relaxation time 
lattice-Boltzmann solverlBSS using the D3Q19 model. The 
velocity field is then used in the advection-diffusion equation 
Eq|3] which is solved by a finite-difference method at the same 
timestep. We assume periodic boundary conditions along the 
flow and vorticity directi ons a nd Lees-Edwards sliding peri- 
odic boundary conditions^^EIl along the direction of the ve- 
locity gradient. FollowingEl we introduce additional internal 
Lees-Edwards planes to overcome an upper limit on shear rate 
that would otherwise be inversely proportional to the system 
size in the gradient direction. (This restriction arises because 
fluid velocities must be small compared to the sound speed in 
the frame of the lattice; the additional planes create a stack 
of lattice slabs whose velocity is nearly matched to the local 



fluid velocity.) The total number of Lees-Edwards planes was 
varied according to the size of the simulation so that the sep- 
aration of two planes was always 16 lattice sites. We relegate 
further comments on the implementation of the Lees-Edwards 
boundary conditions to an Appendix. 

Structural information in Fourier space was obtained from 
the static structure factor C(k), defined as 

C(k)=^(k)^(-k), (8) 

with (|)(k) the normalized Fourier transform of the order pa- 
rameter. The transform itself was calculated using the FFTW- 
library'*^. In principle an ensemble-average or time-average 
should appear in Eq.[8]so that C becomes a smooth function of 
wavevector; in practice an instantaneous measurement, plotted 
as a colour intensity map, was found sufficient. 

For the analysis of the relation between rheological quanti- 
ties and the morphology of the system we need to quantify the 
number of defects. For this purpose we defined an alignment 
function X{r), which measures the average orientation of the 
order parameter gradients, as follows 

^ ^ V'ir)Jv'ir) ^|V^(r)||V^(r')|; ' 

Here V' is a small volume that contains r; in practice V' can 
be chosen to contain r and its nearest neighbours only. 

For ordered regions with small curvature of the lamellar in- 
terfaces, X is very close to unity. Wherever lamellae termi- 
nate at edge dislocations, and in regions where their orienta- 
tion changes abruptly, X departs from unity, taking smaller 
positive values. By setting a threshold X = ^ it is possible 
to detect regions where the order parameter structure deviates 
significantly from what one would recognize visually as an 
aligned lamellar configuration. We found a threshold value of 
t, = 0.8 to be well suited to our requirements. This protocol 
does not directly count the number of defects, since at each 
defect more than one site has an alignment function below the 
threshold. However, as the number of sub-threshold sites per 
defect is reasonably constant, the degree of order under flow 
can usefully be quantified by defining a "defect density" as the 
system-wide fraction of sub-threshold sites: 

pz)(r) = i|rfr'e(^-X(r')). (10) 

Here 9 is the Heaviside function. The conversion factor be- 
tween po and the defect number density (in 2D) or line den- 
sity (in 3D) depends on the lamellar wavenumber q and the 
choice of threshold parameter ^; we do not make such conver- 
sions here and consider po itself as an appropriate measure of 
disorder. 
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3 Parameter Mapping to Physical Units 

The Landau-Brazovskii model belongs to the class of phe- 
nomenological models of liquid-crystalline materials. This 
means it features an isotropic-to-lamellar phase transition 
without explaining the underlying microscopic mechanism. 
Hence there is a many-to-one mapping between physical pa- 
rameters and those of the phenomenological model. However, 
in the context of any particular system, it is possible to indi- 
cate physical parameters for which the chosen model is a good 
representation, and we do this below. 

The parameter space of the model is quite large. Although 
in exploratory work we performed studies across a wider range 
of parameter values, for simplicity in the results presented here 
we restrict attention to fixed values of K = —6 x lO^^LBU and 
c = 7.6 X lO^^LBU (here LRU denotes LB units as defined 
below). These are chosen to give a lamellar wavelength of 
L = 10 LRU (that is, q = ^|k|/2c ^ 27i/10LBU) which offers 
suitable spatial resolution. Moreover the microphase separa- 
tion point, set by flt /Ac — 1.2 x lO^^^LBU, is then a fixed 
constant. 

As usual we set the lattice spacing /, time step Af and mass 
density of the fluid p all to unity to define LB units (LBU). We 
need to calibrate these units (or equivalently, length, energy 
density and time) against real physical units. We do this first 
for thermotropic smectics, then consider the lyotropic case. 

To calibrate the LB length unit, we equate L the lamellar 
wavelength, with that of a typical thermotropic smectic: 10 
LBU — L — 4nm. Hence the LB length unit is about 0.4nm; 
and lm= 2.5 x lO^LBU. Free energy density parameters such 
as Uc — lO^^LBU should then be of order a^kgT /L^ where 
a = 10 is a molecular aspect ratio. Thus the LB unit of energy 
density is about 6 x lO^^Pa, i.e. 1 Pa= L6 x 10"" LBU. 

Turning now to the fluid viscosity r|, after pilot runs span- 
ning a range of several decades, we selected 8.33 x lO^^LBU 
for the production runs reported here. This represents the 
baseline viscosity of the isotropic phase close to the transi- 
tion (not the higher, effective viscosity of the smectic once 
formed). Requiring this to be about ten times the viscosity of 
water (which is J\w = 8.9 x lO^'^Pa s) we find that the LBU 
unit of viscosity is about 0. IPa s. Given the previous result for 
the energy density, we infer that ls= 6 x lO" LBU. The LBU 
timestep is thus around 1.6 x 10" '^s. 

As detailed previously, an estimate of the molecular diffu- 
sivity is D = 2Mb; we set M ~ 0.25LBU thi-oughout this study, 
so that D Ues in the range 2.5 x 10"^ - 2.5 x lO^^^LBU. Tak- 
ing a mid-range value of 5 = 10"'*LBU and using the above 
conversion factors (1 m= 2.5 x 10''LBU, Is = 6 x 10"LBU) 
givesD~ lO^'^m^s"^ which is a reasonable mesogenic value. 

Recall that the mass density of our fluid is unity in LBU. 
However, combining the above conversion factors for length, 
energy density, and time shows the physical and simulation 



units of mass density to be related via 1 kg m ^ = 1 Pa s^ 

m" 

-2 ^ lo-^LBU. Thus our simulations address a fluid of 
mass density lO^'kg m"^ which is 10^ times larger than the 
true value for a typical material. The reason for this is that 
the LB algorithm uses fluid inertia to update the dynamics 
and the discretization should therefore involve replacing the 
real fluid with the densest one possible, while ensuring that 
the Reynolds number Re ^ pyL^/rj (which for real smectics 
is exceedingly small) remains small enough for the simula- 
tion to faithfully describe the low Re limit. As discussed else- 
where ''^^^ Re values up to about 0.1 are then admissable; our 
simulations comply with this requirement. 

Finally, our shear rates lie in the range 10"^ - lO-'^LBU, so 
that the (dimensionless) Peclet number, Pe = yL^/D, which 
expresses the relative importance of advective to diffusive 
transport, lies in the range 5 — 800. These are moderate to 
high values where significant structural rearrangement can be 
expected. However, for thermotropic smectics, given our time 
unit of 1LBU= 1.6 x lO^^^s, they correspond to strain rates 
of order 10''s"^ This lies well above the range typically 
addressable with rheological experiments, but below that re- 
cently given in a theoretical study 

The shear rates needed to achieve this range of Pe fall dra- 
matically when lyotropic smectics are considered. For ly- 
otropics one has L ^ 40nm; the LB length unit is thus about 
4nm. The factor a in the free energy density is now of order 
unity (since the main contributions are entropic and for flexi- 
ble membranes there is about one relevant degree of freedom 
per cube of side Lj^. The LB unit of energy density is then 
about 6 X lO^Pa, i.e. 1 Pa= 1.6 x 10"*' LBU. Both the vis- 
cosity r| and the molecular diffusivity D remain comparable 
(within one order of magnitude) to the thermotropic values; 
from the viscosity we then determine that ls= 6 x lO^LBU or 
ILB timestep = 1.6 x 10"^s. Accordingly the strain rates rel- 
evant to the Pe range addressed here are of order 10^ s"' and 
well within experimental reach, but below that addressed in 
molecular dynamics simulations'^S^. Note that for lyotropics 
we are effectively simulating a fluid of even higher mass den- 
sity than for thermotropics; but since the simulated Reynolds 
numbers are unaffected by the change in parameter mapping, 
they are still sufficiently small"*^. 

To further simplify our parameter exploration, throughout 
what follows we restrict attention to states with b = —a. This 
convention has been widely used to study bulk spinodal de- 
composition in binary fluids where it can be absorbed in a 
rescaling of the order parametei^l put differently, b = —a 
ensures that coexisting bulk ordered phases have (|) = ±(|)o 
with (|)o = 1. In lamellar states of b = —a the amplitude 
of the lamellar ordering need not be unity but varies like 
(|)^ ^ -y/l +ac/b as stated previously. Thus the chosen value 
of b can be thought of as controlling whether the peak order 
parameter variation is close to that arising for bulk coexistence 
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(when b ^ a^-), or whether the amplitude (|)^ is instead much 
larger (b <^ Oc) and determined by a competition between the 
gradient terms and a relatively weak quartic nonlinearity. In 
this study we choose three values, b = 5x 10^"*, 1 x lO^'* and 
5 X lO^^LBU so that the dimensionless quantity (3 = b/ac is 
respectively 4.2,0.83,0.42, lying near the crossover between 
these regimes. The observed (|)^ range is about 1.2— 1.5, 
somewhat narrower than the range 1.1 — 1.8 predicted from 
(^q ^ ^J\A-~a^Jb. However the latter estimate, which does not 
treat the nonlinear term exactly, is only a guideline. 

Alongside the Reynolds number (which we deem irrele- 
vant), P, and the Peclet number, we can identify two further 
dimensionless numbers of interest. One is the Schmidt num- 
ber Sc = V[/pD, which is large when momentum can diffuse 
much more rapidly than molecules - as applies in almost all 
dense fluids. In our simulations Sc is of order 100 — 1000 
which is safely large. (The work of KumaranIS! instead ad- 
dresses Schmidt numbers in the range 0.2 — 5 where differ- 
ent physics may apply.) The second parameter controls the 
strength of feedback between distortions of the (|) field and the 
fluid motion. This can be expressed in terms of a characteris- 
tic value of the Peclet number, Pe*, at which the free energy 
density in the order parameter field that couples to deforma- 
tion (Eq|2|l, g = (qc + b)ac/4-b, is comparable to the viscous 
stress Oy = iiy = {r\D/L?')Pe. Equating these gives 



Pe* 



{ac + b)ac 
8P rjM' 



(11) 



With our chosen parameters (in LBU, = 1.2 x 10 = 
10,r| = 8.3 X 10-2, M = 0.25) we find Pe* - 180 at p = 4.2 
and Pe* - 5000 at P = 0.42. The maximum values of Pe/Pe* 
that we address below are of order 0.2 and 0.08 respectively 
in the two cases. 

For later convenience we may define a new dimensionless 
number, which we can call a capillary number. (The capil- 
lary number is conventionally defined as the ratio of viscous 
to interfacial forces in binary fluid flow.) We define this as; 



Ca = 



Pe 

Pe* 



4r[jb 



{ac + b)ac 



(12) 



We will argue below that this could be a more relevant choice 
than Pe for understanding the dependence of morphology on 
flow rate in smectic systems, at least when Pe* ^ 1. 



4 Results and Discussion 

We first present detailed results for two parameter sets (P = 
4.2,0.42) in two dimensions which we call systems A-2D and 
B-2D. System lattice sizes and shear rates are reported in Table 
[T| all remaining parameters are the same in both systems as 
detailed above. A third parameter set (p = 0.83, system C- 
2D) gave results broadly similar to B-2D. We will later report 
simulations in three dimensions for systems A-3D and C-3D 
which share parameters with A-2D and C-2D (see Table [T]). 
Note that our simulations in two dimensions are equivalent 
to assuming, in three dimensions, that the system maintains 
translational invariance in the vorticity direction. We find this 
is an increasingly good approximation at very low shear rates. 

4.1 Two Dimensions 

Before presenting specific results for the parameter sets de- 
tailed above, we briefly describe the effect of wider param- 
eter variations. Within the Landau-Brazovskii free energy 
functional (Eq. [TJ a range of parameter variations can pro- 
duce stable lamellar phases, as explored in previous stud- 
jgg ,25,28|4i,52j53, , 'pj^g nonequilibrium morphology of the lamel- 
lar structure, and the phase ordering kinetics (starting from 
a slightly noisy, uniform initial state) do however depend on 
the choice made for these parameters, and also for the viscos- 
ity. Formation and ripening of the lamellar structure is slower 
when P is small, consistent with the reduced initial driving 
force for ordering in this case. Higher viscosities lead to more 
locally disordered structures in which there is a shorter sepa- 
ration between kinks or edge defects as one moves along the 
lamellae; conversely lower viscosities typically yield "longer" 
lamellae. This effect has been reported before we found it 
to be more pronounced for softer (lower p) systems. At fixed y, 
higher values of r| were found to cause greater inhomogene- 
ity of lamellar shape and thickness, consistent with the fact 
that Caocr[ (the larger the shear stress, the larger the lamellar 
deformation). 

We turn now to system A-2D (p = 4.2). The top panel 
in Fig. [T] shows the order parameter pattern (|) at the end of 
the equilibration phase prior to shearing. Although the true 
equilibrium structure would show homogeneous alignment, 
this cannot be reached in a reasonable time scale unless shear 
is applied. The resulting metastable structure is amorphous 
and isotropic, but locally exhibits the same lamellar spacing 
L = 10 LBU as the fully ordered equilibrium state. If shear 
flow is imposed, the lamellae start to align in the flow direc- 
tion, which induces a morphological changes to the structure. 
A typical snapshot in the resulting steady state at Pe = 40 
{Ca = 0.22) is shown in the middle panel of Fig. [T] In this and 
all subsequent images, the applied flow velocity is horizontal 
and the velocity gradient vertical, with the flow directed to the 
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ran # 




time /lO^ 




Ca 


Pe 


lO^Y 


A 














1 


1024 X 512 X 1 


0-3.2 


50 











2 


1024 X 512 X 1 


3.2-6.4 


50 


0.224 


40 


10 


3 


1024 X 512 X 1 


6.4-9.6 


50 











4 


512 X 256 X 1 


0-3.2 


50 











5 


512 X 256 X 1 


3.2-6.4 


50 


0.448 


80 


20 


6 


512 X 256 X 1 


3.2-25.6 


50 


0.224 


40 


10 


7 


512 X 256 X 1 


12.8-38.4 


50 


0.112 


20 


5 


8 


512 X 256 X 1 


51.2-76.8 


50 


0.056 


10 


2.5 


9 


512 X 256 X 1 


51.2-102.4 


50 


0.028 


5 


1.25 


10 


256 X 128 X 1 


3.2-16.0 


50 


0.22 


40 


10 


11 


256 X 128 X 128 


0-6.4 


50 











12 


256 X 128 X 128 


3.2-6.4 


50 


0.448 


80 


20 


13 


256 X 128 X 128 


3.2-9.6 


50 


0.224 


40 


10 


14 


256 X 128 X 128 


6.4-12.8 


50 


0.112 


20 


5 


B 














15 


1024 X 512 X 1 


0-3.2 


5 











16 


1024 X 512 X 1 


3.2-9.6 


5 


0.082 


400 


10 


17 


1024 X 512 X 1 


9.6-12.8 


5 











18 


512 X 256 X 1 


0-3.2 


5 










19 


512 X 256 X 1 


3.2-9.6 


5 


0.164 


800 


20 


20 


512 X 256 X 1 


3.2-25.6 


5 


0.082 


400 


10 


21 


512 X 256 X 1 


12.8-38.4 


5 


0.041 


200 


5 


22 


512 X 256 X 1 


25.6-64.0 


5 


0.020 


100 


2.5 


23 


512 X 256 X 1 


51.2-76.8 


5 


0.010 


50 


1.25 


C 














24 


256 X 128 X 1 


0-3.2 


10 











25 


256 X 128 X 1 


3.2-9.6 


10 


0.252 


400 


20 


26 


256 X 128 X 1 


3.2-12.8 


10 


0.126 


200 


10 


27 


256 X 128 X 128 


0-3.2 


10 











28 


256 X 128 X 128 


3.2-9.6 


10 


0.252 


400 


20 


29 


256 X 128 X 128 


3.2-12.8 


10 


0.126 


200 


10 



Table 1 Simulation parameters: After initialization from a noisy 
uniform state, the systems were first equilibrated at zero shear rate, 
followed by a sequence of various shear rates, ranging from 
Y= 1.25 X 10-^ to 2 X lO-^LBU. In some cases we performed a 
final switch off run at zero shear rate. The other model parameters 



IX 10^ 



: 7.6 X lO-'^.M = 0.25 and r| = 8.33 x 10 



-2 



were the same throughout all simulations. 



right in the upper half and to the left in the lower half of the 
picture. After a short transient during which the order param- 
eter pattern remains relatively homogenous (not shown) there 
emerge extended regions in which the smectic layers are either 
expanded or compressed, so that L locally deviates from the 
equilibrium value. These mesoscopic inhomogeneities are in- 
£lined at an angle to the flow and remain continuously present 



J£_Lin the steady state, but are not stationary: their local features 
_are advected by the flow and evolve continuously by deforma- 



tion. 




Fig. 1 Order parameter pattern (|)(r) of system A-2D (thresholded to 
enhance visibility): after equilibration (timestep / = 3.2 x 10^, top), 
in steady shear at Y = 1 x 10^'* (f = 5 x 10^, centre) and after 
switch-off (/ = 1 .28 X 10^, bottom). Gray scaling from black to 
white corresponds to values of — 1 < (]) < 1 . 



We found that after the shear flow was switched off the 
lamellar spacing returned to its equilibrium value everywhere. 
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with expanded regions contracting and vice versa. The final 
state, after all relaxation has (virtually) come to a standstill, 
is shown in the bottom panel of Fig. [T] Most of the lamel- 
lae remain oriented along the flow direction, but depending 
on the local expansion and compression undergone during the 
previous flow, some regions retain strong misalignment which 
does not anneal away. The separation between defects as one 
moves along the lamellae varies considerably in space; both 
short and long separations can be seen. This applies both to 
the aligned layers and to layers that are misahgned with the 
flow. 




Fig. 2 Order parameter (])(r) of system B-2D: after equilibration 

(timestep t = 3.2 x 10^, top), in steady shear at Y = 1 x 10^^ 

(f = 8 X 10^, centre) and after switch-off (t = 1.28 x 10^, bottom). 

Fig. [2]shows equivalent states for system B-2D (p = 0.42). 
The Peclet number during the shearing phase is now Pe = 400 
(Cfl = 0.082). All the generic results of Fig. [T] are the same. 



but the defect density, both during flow and after switchoff, is 
significantly reduced. (However the bands of misalignment in 
the flowing state are, if anything, even more apparent.) Conse- 
quently, the lamellar alignment along the flow direction after 
switch off is also more complete. 

It is notable that at the shear rate shown in these figures, the 
more ordered system is B-2D which has higher Pe than A-2D 
but smaller Ca. Varying y within system A-2D or B-2D es- 
tablishes unambiguously that the steady state order decreases 
with shear rate at fixed thermodynamic parameters (Fig|6]be- 
low). This points strongly to Ca as the more relevant measure 
of flow rate when discussing onset of the disordered phase. 

In principle the dependence of the state of order on flow 
rate could be smooth, or could entail a nonequilibrium phase 
transition - for instance if the perfectl y ordered state ceased 
to be stable above a certain shear rateSEHlIllll in practice, 
however, starting as we do from an initially amorphous lamel- 
lar texture, a perfectly ordered state is generally not reached 
for large systems even if it might be stable once formed. In 
any case it seems likely that at high Pe, any defects introduced 
into the smectic state by deformation of an initially imperfect 
structure (whose layers are then dilated and compressed by 
the flow) cannot anneal away or annihilate each other rapidly. 
Such irregularities then build up until the steady state defect 
density is determined by a balance between their rate of intro- 
duction and rate of annealing. 

Recently Kumaran and Raman studied a single shear rate 
and varied the Schmidt number Sc and system size, reporting 
enhanced alignment on decreasing Sc and in small systems'^. 
Our comparative studies on systems A and C have confirmed 
both trends (although, as mentioned previously, our Sc values 
are much higher). The same authors also discuss a distinct 
mechanism of defect creation caused by layer compression 
and splitting. Although we have not checked this in detail, 
this may be consistent with our observation of a strongly com- 
pressed and heterogeneous morphology, with bands oriented 
normal to the compression axis, as shown in Figs. [T] and |2] 

A somewhat different morphology in sheared lamellar sys- 
tems was reported by Xu et al.^**. This morphology exhibits 
macroscopic inhomogeneity between the centre of the simu- 
lated domain and the regions near the walls, accompanied by 
variations in the local shear rate. This contrasts with our own 
work which shows translationally invariant statistics in all di- 
rections. These differences may be a result of Xu et al.'^Susing 
mixed boundary conditions, as discussed already in the intro- 
duction, in which hard walls for the fluid flow are combined 
with periodic BCs on the order parameter field. 

To quantify our observed morphology more closely, we use 
the alignment function X introduced in Eq. [9] Fig. [3] shows 
snapshots of the time evolution of system A-2D at modest 
shear rate. All lattice sites with X{r) < 0.8 are marked in red 
and they obviously coincide with lamellar endpoints (edge dis- 
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Fig. 3 Defect stmcture of system A-2D: at timestep ? = 5.12 x 10^ 
(top) and 1.024 x 10^ (bottom) with an imposed shear rate 
Y= 1.25 x 10~^. The system has been pre-sheared at stepwise 
decreasing rates (see Table[TJ. Sites where the alignment function X 
is smaller than 0.8 are highlighted red. The system size was 
WvXA'yXA', = 512x 256 X 1. 



locations) or other local defects. The upper panel gives the sit- 
uation at the lowest applied shear rate y = 1 .25 x 10^^, imme- 
diately following a pre-shearing protocol in which the strain 
rate was ramped down stepwise from higher values (see Table 
[T]l. While large aligned regions are present there are other re- 
gions of high defect density. The lower panel of Fig.|3]shows 
a later state with a total strain incremented by a further 6400% 
lattice units (i.e. by a displacement 64 times the system size). 
Obviously many defects have annealed away and the smectic 
layers are now rather well aligned with the flow direction. 

Defect annihilation to this large extent is only found at 
low enough shear rates. A prerequisite seems to be that 
the defects can smoothly approach each other and then an- 
nihilate. This is prevented by the buckling of layers above 
a certain shear rate; the undulations so formed ultimately 
pinch off and break, resulting in the creation of new defects. 
The onset of defect proliferation arises at characteristic flow 
rates % = 1 .25 — 2.5 x lO^^LBU in system A-2D and around 
Yc = 1 - 2 X lO^'^LBU in B-2D. (Here the subscript c could 
stand for "critical" or for "crossover" depending on whether 
there is a singular transition or not; see above.) These values 
are for systems of size 512 x 256 x 1; we found evidence of 
an increasing trend with system size, but the extremely long 
run-times involved (^ 10^ LB timesteps) precluded quantita- 
tive exploration. These results correspond respectively to val- 
ues of the Peclet number, Pe^ obeying 5 < Pec < 10 (system 
A) and 400 < Pec < 800 (system B). However, in terms of 
Ca = Pe/Pe* one has 0.056 < Cac < 0.112 (system A) and 
0.082 < Cac < 0. 164 (system B). Given the 80-fold difference 
in Pbc, the similarity of Cac in the two different systems is 
remarkable. 

Boldly extrapolating across parameter regimes that we have 
not explored in any detail, our results motivate the hypothesis 
that the onset of defect proliferation occurs generically in 2D 
at Cac ~ 0.08. Note that Ca, which comprises the ratio of a 
viscous stress to an elastic modulus, can also be viewed as a 
strain. In that context one can imagine that defects prolifer- 
ate when the steady-state elastic deformation of the smectic 
exceeds a yield strain value (for which a threshold of a few 
percent is reasonable). This makes sense if, at Ca ~ Cac, the 
stress contributed from the order parameter is (a) comparable 
to the viscous stress from the background fluid - essentially 
the force balance argument made above - and (b) primarily 
elastic in character. The latter in turn requires that Pec is large 
and hence that Pe* is larger still. A different mechanism might 
apply at small Pe*, since in that case layers relax rapidly by 
diffusion on the time scale set by Cac- 

In Fig. |4] we explore the correlation between the inhomo- 
geneous morphology in the flowing state and locally defined 
rheological quantities. The pictures show the order param- 
eter (|)(r), shear stress density Oxy{r) and first normal stress 
difference density Ni{r) — O^^ii') ~c^>t('") of system A-2D 
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Fig. 4 Order parameter (|)(r) (top), shear stress CTtv(r) (centre) and 
first normal stress differences Ni (r) (bottom) of system A-2D at 
f = 4 X 10^ (see TablefTl. 



in steady shear flow at shear rate y = 1 x 10^^ (Pe = 40, 
Ca = 0.22). Here and below we adopt a coordinate system 
where x is along the flow direction and y the velocity gradient. 
In regions where the lamellar width is larger than in the equi- 
librium state, shear stresses and first normal stress differences 
show low positive or negative values. Large positive values of 
both quantities coincide with regions where the lamellae are 
under compression. Peak values (in LBU) are in the region of 
—2 X lO"'* and 3 x 10""* for shear stress and —3.5 x 10"^ and 
7 X 10^"* for the normal stress differences. Despite having lo- 
cally negative regions, spatial and ensemble averages of these 
quantities are always positive. Very similar results, but with 
smaller peak magnitudes, are obtained for system B-2D at the 
same shear rate. 

Fig. |5] shows structure factors C{k) for the same states as 
Fig. [T] C{k) is initially isotropic after layer formation and 
equilibration in the system at rest with, as expected, a distinct 
peak at |*| ~ ±271/10 LBU (recafl L = 10 LBU). In steady 
shear flow the structure factor develops two distinct but rather 
broad peaks at fc^ ~ ^0.0l2K,ky ~ ±0.227t. Thus the inho- 
mogenous morphology of the lamellar structure is clearly en- 
coded in the structure factor. The wavevector at the peak has 
a small component in the x-direction, indicating the average 
tilt of the lamellae away from the flow direction towards the 
elongational axis of the flow (which in shear is at 7t/4 to the 
flow direction). 

When the shear flow is switched off, the system relaxes 
again to its equilibrium layer spacing and the small tilt of the 
C{k) signal decreases to an undetectable level: this is consis- 
tent with the fact that, in Figs. [T] and |2] (bottom panels) it is 
easy to see when the system has been sheared horizontally, 
but not in which sense. The intensity of the peak is slightly 
enhanced in the final state due to the improved orientation of 
the layers along the flow direction. 

Fig. |6] shows the time evolution of the defects density po 
of system A-2D in the sequence of stepwise decreasing shear 
rates alluded to previously (starting at Y= 2 x 10^"* and rang- 
ing down to Y = 2.5 x 10^^, Table [lb. The runs with the 
two highest shear rates were started from the amorphous, 
post-equilibration configuration, whereas all others were pre- 
sheared using the final configuration of the preceding run as 
the initial state. The total strain between the beginning and 
end of every individual run was at least 3 x 10^%. Interest- 
ingly, following each step down in shear rate, both systems 
A-2D and B-2D continuously anneal to lower defect densi- 
ties on a timescale ~ lO^LBU. If the shear rate is turned up 
again (not shown), the lamellae buckle and new defects re- 
cur on a timescale corresponding to a stain of about 100%. 
This finding offers further support for a dynamical equilib- 
rium between defect creation and annihilation above a critical 
(or crossover) shear rate, below which defects progressively 
annihiliate and near-perfect alignment eventually results. To 
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Fig. 5 Structure factor C(k) of system A-2D: Shown are states after 
equilibration at f = 3.2 x 10^ (top); in steady shear at y = 1 x 10^"* 
and f = 5 X 10^ (centre) and after switch-off at f = 1.28 x 10* 
(bottom). The sections shown correspond to wave vectors on the 
interval kf,ky G [— 7t/2,3t/2] LBU with A;, and ky along the 
horizontal and vertical coordinate direction, respectively. 



test this assumption we performed simulations starting from a 
very well aligned state and increased the shear rate. The defect 
density indeed goes up (not shown) and reaches similar values 
to before. Note however that the defect annealing time (of or- 
der lO^LBU) is quite long compared to any inherent relaxation 
time; also, the final state after switchoff strongly depends on 
the previous flow. Thus the rheology and microstructure of 
these materials is inherently history dependent. 
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Fig. 6 Defect density po vs time: The main picture shows data for 
system A-2D in a sequence of simulations with consecutively 
lowered shear rate y = 2 x 10^^ (red solid), 1 x lO^'* (green 
long-dashed), 5 x 10^-* (blue short-dashed), 2.5 x 10^^ (magenta 
dotted) and 1.25 x 10^^ (cyan dashed-dotted). Horizontal lines 
indicate time averages after entering the apparent steady state 
(where present). The inset gives data for the system A-3D and shear 
rates Y = 2 x 10""*, 1 x 10""* and 5 x 10"^ (left to right). 

The top picture in Fig. |7] shows the spatially averaged total 
shear stress (o.q (r)) = Oxy At low shear rates, stress fluc- 
tuations have longer temporal correlations but smaller ampli- 
tudes. (The latter dependence is obscured by the logarithmic 
vertical scale in this plot.) The root-mean-square temporal 
fluctuations can be reduced by increasing system size (by a 
factor 3 or so for the runs shown in Fig. [T] and Fig. |2|. In 
contrast, the averaged shear stress is almost independent of 
the system size, which makes it possible to simulate longer 
times using smaller system sizes while still giving accurate 
predictions for the bulk rheological behaviour. We pursue this 
strategy here. 

System A-2D shows clear shear thinning behaviour; val- 
ues of its apparent viscosity at various shear rates are given 
in the caption of Fig. 7. The corresponding range of Ca is 
0.028-0.448. Using Eq|5] we can identify the stress contri- 
bution arising from the deviatoric part of the order parameter 
pressure tensor {—P„{r)) = — P,j.; this is shown in the centre 
graph in Fig. |7] The temporal fluctuations are in phase with 
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Fig. 7 Top: average total shear stress a^y versus time for system 
A-2D: The horizontal lines are steady state time averages. For shear 
rate Y= 1.25 x 10^^2.5 x 10"^5 x 10"^ 1 x 10""*, 2 x 10""^, we 
find r\app/^ = o.vv/t|Y= 2.49,2.25,2.18, 1.93, 1.74; centre: average 
order parameter contribution to the stress tensor — Piv vs time for 
system A-2D. The inset shows data for system A-3D and shear rates 
Y = 2 X 10""*, 1 x 10"^ and 5 x 10"^ (left to right); bottom: average 
chemical contribution to the pressure tensor —Pxy vs. defect density 
for system A-2D. 



those of the total stress and decrease with decreasing shear 
rate. The bottom image in Fig. |7]plots —Pxy against the de- 
fect density p^; this shows a linear relation across almost the 
whole range of shear rates. 

In contrast, system B-2D does not exhibit clear 
shear thinning; the apparent viscosities obey ^app/^ — 
1.57,1.62,1.55,1.42,1.87 for the same set of shear rates as 
used in Fig.|7] To within the measured accuracy, system B-2D 
is thus rheologically quasi-Newtonian. However, the range of 
Ca in this case is lower, 0.01-0.164, and the order parameter 
stress creates only a relatively modest uplift to the baseline 
viscosity. More significantly, Ca < Ca^ for all but the highest 
applied shear rate so that the steady states found in system 
B-2D are far more ordered than in A-2D. This results from a 
defect annihilation process, occurring in the early stages of 
ramp-down, which we found to be accompanied by a sharp 
drop in P„ (not shown). 
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Fig. 8 Spatially averaged first normal stress differences 
Ni = Oxx — Oyy vs time for system A-2D. The inset shows data for 
system A-3D and shear rates Y = 2 x 10^"*, 1 x 10^^ and 5 x 10^^ 
(left to right). Horizontal solid lines are time averages. 

Fig. |8] shows averaged first normal stress differences 
(A^i(r)) = A^i versus time for system A-2D, where we report 
the same stepwise decreasing flow protocol as before. For the 
first several steps A^i decreases monotonously with decreasing 
shear rate. However, at the lowest shear rates the behaviour is 
not monotic: A^i eventually rises again. (The resulting time av- 
erages, shown as horizontal lines, do not therefore reflect true 
steady states.) This somewhat surprising feature has a simple 
explanation. The time sequence of the defect density in Fig. |6] 
reveals that the rise in A^i is correlated with a drop in po, and 
linked to the onset of defect annihilation (which happens only 
at low shear rates). When two defects meet and cancel, the 
local lamellar width alters by a certain amount. This causes 
frustration as the system is unable to snap back locally towards 
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its equilibrium layer spacing. Hence the normal stress differ- 
ence Ni — Gxjc — Oyy fails to decay to zero, and remains pos- 
itive at large times, as appropriate for a system of metastable 
and (on average) compressed smectic layers. The time depen- 
dence of the normal stress difference is also anomalously slow, 
whereas the shear stress reaches their steady state values more 
quickly. The relatively large steady state value of A^i and this 
slow dynamics may partly be a finite size effect caused by the 
relatively modest number of smectic layers within the simu- 
lated volume. The normal stress in system B-2D (not shown) 
behaves differently: in this case A^i falls monotonically during 
the defect annihilation phase. It is not clear whether this differ- 
ence should be attributed to the change in P,Ca or Pe between 
the two systems, or indeed a combination of those changes. 
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Fig. 9 Flow curves: Time-averaged shear stress (cT.iy) and first 
normal stress differences (A'l ) vs shear rate y for system A-2D and 
B-2D. The dash-dotted lines are power-law fits to the data. 

Fig. |9] shows 'flow curves' representing the shear rate de- 
pendence of time averaged stresses o„ and A^i . This is cal- 
culated from the stepwise ramp data as detailed previously, 
so the A^i results at low shear rates should not be taken as 
definitive. System A-2D shows shear thinning as already dis- 
cussed, apparently following a power law a^y = with ex- 
ponent m ~ 0.87. A recent experimental and theoretical work 
studies the nonlinear relation between shear stress and shear 
rate in lyotropic lamellar phasesES interestingly the authors 
identify the motion of screw defects as being responsible for 
the layer tilting, folding and formation of defects. They re- 
port an experimental value of m ~ 1/1 .44 ~ 0.694, which is in 
good agreement with their theoretical prediction of m — 2/3. 
Our value of m ~ 0.87 for the disordered flow regime is larger. 
One possible explanation for this discrepancy is that screw de- 
fects, which do not exist in our 2D simulations, play a key role 
for this theory. On the other hand our fitted power law must 
be treated with great caution since the fit includes data from 



Y both above and below y^.. It is somewhat surprising that the 
morphological changes around 1.25 x 10^^ < % < 2.5 x 10^^ 
in system A give no visible anomaly or feature in the flow 
curves. For system B-2D the data is consistent with an ex- 
ponent m = 1, indicating more or less constant viscosity and 
no pronounced shear thinning (as previously described). This 
appears to be due to the decreasing number of defects at shear 
rates below y < 2 x 10^^; the order parameter stress P„ is then 
dominated by the background viscosity term r|y. 

In neither system does 0„{j) offer any clear signature of a 
yield stress. There may be two ways to understand this. Ei- 
ther, as detailed previously, this is because at the lowest shear 
rates studied here the selected state is a fully aligned one with 
the layering wave vector perpendicular to the flow: such a state 
should have a purely viscous response. However, a defect tex- 
ture may also lead to a viscous response in the linear regime 
if flow occurs via permeation ' (i.e. by fluid flow not altering 
appreciably the order parameter pattern). 

To further characterise the small shear rate regimes in our 
simulations, we have performed additional simulations on 
fully aligned states, with different layer spacings. While the 
shear stress contribution in this case is negligible as y -> 0, 
A^i does not vanish for vanishing shear rate, but depends sen- 
sitively on the nonequilibrium layer spacing in the range of 
shear rates addressed here. This can be described as a "quasi- 
Newtonian" behaviour. As expected, A^i is positive for over- 
compressed layers, and negative for overdilated layers. The 
values of A^i obtained in system B-2D (see Fig. 9) are compa- 
rable to those observed in a quiescent lamellar smectic system 
prepared by compressing the layers by about 5-10%. 

4.2 Three Dimensions 

The computational cost of repeating all the above studies in 
three dimensional systems of equivalent size remains pro- 
hibitive, particularly when addressing steady state properties 
which emerge only after a million or more LB timesteps. We 
therefore focus in 3D on smaller systems (Table [T]i, but still 
large enough in principle to show qualitatively similar struc- 
ture under flow. Our studies use two sets of parameters. Sys- 
tem A-3D is identical to A-2D except for the introduction of 
the third dimension. System C-3D has parameters interme- 
diate between A-2D and B-2D (Table [TJ and was chosen in 
preference to the latter because B-2D ordered completely if 
the system size was matched to a slice of B-3D. (We see be- 
low however that in fact this is a rather general outcome in 
3D.) 

Fig. 10 shows the post-equilibration amorphous structure 
and the steady-state structure under flow in system C-3D (A- 
3D is similar) at system size A^^ x Ny x A^, = 256 x 128 x 128. 
The latter is compared with C-2D and A-2D in a system of 
the same x and y dimensions (256 x 128 x 1). The quiescent 
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Fig. 10 Order parameter (|)(r) in two and three dimensions: Top pictures sliow system C-3D after equilibration (left) and in shear flow (right, 
Y = 2 X 10^'*) at f = 6.4 X 10^ LBU. The opacity of one phase has been reduced to enable a better view into the lamellar structure. The bottom 
left picture shows system C-2D at the same time step and shear rate. The bottom right picture is system A-2D at y = 1 x 10^^ and time step 
t = 9.6 X 10^ LBU. Results for A-3D are morphologically indistinguishable from those for C-3D. In all these images the flow direction is 
horizontal, to the right in the upper half of the sample and to the left in the lower half. 
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Fig. 11 Structure factor C(k) of quiescent system A-3D, showing a 
cut along kx = after equilibration at f = 3.2 x 10^. Shown are 
wave vectors on the interval ky,k. e [— 7l/2,n/2] with ky and k, 
being the horizontal and vertical direction, respectively. Data for 
C-3D (not shown) is virtually the same. 



metastable state of system C-3D after equilibration is more or 
less the same as it would be in two dimensions. Indeed, a slice 
through the structure factor along kx = (cuts along ky,k^ = 
look identical) is shown in Fig. 1 1 and resembles that in Fig.|5] 
However, under shear at y = 2 x lO""* {Pe = 400, Ca ^ 0.252) 
this structure anneals readily into a regular lamellar structure 
after a transient period of about 6 x 10^ LB timesteps. This 
occurs for even the highest shear rates applied, indicating that 
there is a very strong tendency in three dimensions to form 
ordered lamellar arrays. The same behaviour was seen in sys- 
tem A-3D for which the highest flow rate used corresponds 
to Pe = 80 and Ca = 0.448. (In contrast, the 2D results at the 
same flow rates and system sizes feature the typical fluctuating 
pattern of breaking and merging lamellae that was described 



previously and is visible in the lower panels of Fig. 10 ) The 



structure factor seen in 3D (see Fig. [T2jl is accordingly close 
to that of a perfectly ordered lamellar state, in contrast to the 
corresponding results in 2D (Fig. |5]l. These findings suggest 
that any characteristic strain rate % at which defects prolifer- 
ate must be markedly higher in 3D than in 2D. Specifically, 
if Ca is the relevant control parameter, our results for A-3D 
imply that Ca^ > 0.448 in three dimensions. 

Note however that the order we observe is not completely 
perfect as a mismatch of the layers occurs leading to branching 
points and edges (e.g. at the centre of the ordered 3D-system 
in Fig. [To] ). The Burgers vector is oriented along the gradi- 
ent direction with edges along the flow direction. This forms 
a structure that is translationally invariant along the flow di- 




Fig. 12 Structure factor C(k) of system A-3D in shear flow 
(Y = 1 x 10""*) at t=960k. Results are shown in the planes atkx = 
(top panel), k, = (middle panel), and ky = 37l/16 (middle panel). 
The horizontal and vertical directions are respectively along k, and 
ky (top panel), ky and ky (middle panel), kx and k, (bottom panel). 
The wave-vectors in the horizontal and vertical directions are in 
[— jr/2,7t/2] for all panels. The structure factor is akin to that 
expected of a regular lamellar phase. 
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rection and in the final steady state the location of the edge 
defects does not change. Notably, the corresponding orienta- 
tion of edge dislocations does not exist in two dimensions. 

A closer look at the kinetic evolution of the order parame- 
ter reveals that the additional dimension appears to play a key 
role in promoting the formation of the homogeneously layered 
structure. The ordering process takes place in two consecutive 
steps. Firstly, the equilibrated state undergoes a transition to a 
semi-ordered state, which consists of small patches of layered 
regions, with the layer normals along gradient direction. The 
length of the lamellae is initially quite short, but gets longer 
along the flow direction as the simulation proceeds and the 
structure is advected. It remains, however, relatively short 
along the vorticity direction, with the individual layers being 
disconnected. In a second step the layers merge by linking up 
along the vorticity direction. This is possible as the layers can 
be locally translated along the flow gradient direction, so that 
the average order parameter within a plane can change and is 
not conserved as in the two-dimensional case. 

We next quantify the time system A-3D takes to anneal into 
its final ordered structure starting from an isotropic initial con- 
figuration. (These and subsequent measurements are reported 
for A-3D only; C-3D is qualitatively similar.) The inset of 
Fig. |6] gives the defect density as defined by Eq. 10 versus 
time for shear rates y = 5 x 10^^ — 2 x 10^^. Within a few 
hundred thousand LB timesteps the defect density decreases 
by two orders of magnitude, resembling the behaviour seen 
in system A-2D at shear rates one decade smaller The total 
strain where the curves level off, and ordering is achieved is 
about y= 160 (atY=2x 10"'*),Y= 120 (aty = 1 x lO"'^) and 
Y=65 (atY=5 X 10"^). 

The reduced number of defects is reflected in the magni- 
tude of the order parameter stress contribution the 3D be- 
haviour is shown in the inset of Fig. [8] The final contribution 
to the total stress is created by a small number of remaining 
defects (hence the inconsistent trend of the tiny plateau value 
with strain rate, which is not a concern here) and is quite neg- 
ligible compared to the background viscosity contribution r|y. 
The drop in P„ by three orders of magnitude during the defect 
annihilation process can formally be viewed as shear thinning 
of a quite extreme degree - albeit one masked, in any macro- 
scopic measurement of Gxy, by the Newtonian background. 

First normal stress differences A^i = {Oxx — <^yy) for system 
A-3D are shown in the inset of Fig. [8] For all shear rates ad- 
dressed here system A-3D attains similar values of the nor- 
mal stress; these are approached monotonically and about one 
decade smaller than those seen in system A-2D. The mecha- 
nism reported previously for the buildup of normal stress on 
approaching the ordered state in 2D systems seemingly does 
not apply in 3D, perhaps because of the much more efficient 
pathways to relaxation and mutual annihilation of defects. 



5 Conclusions 

Many earlier studies have proposed a picture in which a 
layer undulation instability in three-dimensional smectic lich 
uid crystals emerges above a critical shear rate ^SiSl 
Molecular dynamics simulations on lyotropic smectics place 
this in a regime of flow r ates w ell above those we have ap- 
plied in the present studyGSUl The critical shear rate for 
thermotropic smectics that has been given in a theoretical ap- 
proach is slightly above ours. The status of this instability 
and its effects on rheology in experiments remains somewhat 
unclear. Also unclear is its status in two dimensions. 

In the present work we have studied the nonlinear rheology 
of lamellar systems in two and three dimensions by means 
of large-scale simulations of a well-established model. Af- 
ter quenching from a homogeneous phase a metastable amor- 
phous lamellar structure appears, similar to previous simula- 
tions in two dimensions using comparable models and param- 
etersl221t2Sl. Also in accord with some previous studies, under 
flow in two dimensions we found the morphology of the order 
parameter to show two distinct regimes: one with very few de- 
fects which arises at low shear rates, and a defect rich phase 
at higher ones. We found that the critical (or crossover) strain 
rate % separating these phases depends somewhat on system 
size, but much more strongly on thermodynamic and kinetic 
parameters. We have argued that this dependence might be 



summarized in terms of a capillary number Ca (Eq. 12 1 that 



is the ratio of viscous stress to the gradient free energy den- 
sity of the lamellar phase. This controls a steady-state elastic 
strain in the smectic which, when it exceeds some characteris- 
tic value Cuc, might cause yielding of layers and defect prolif- 
eration; our estimate is 0.056 < Coc < 0.112 (system A) and 
0.082 < Coc < 0.164 (system B). However we emphasize that 
the evidence for Ca as a control parameter is currently very 
limited and it would take a painstaking exploration of param- 
eter space to confirm this picture. 

A notable feature of the disordered flow phase is our obser- 
vation of striated heterogeneity in the state of smectic order, 
with strongly dilated and compressed regions oriented along 
the extensional axis of the flow. These are present in steady 
state but are nonstationary. The extended regions with lamel- 
lar spacing larger than at equilibrium exhibit either negative 
or small positive shear stress and first normal stress differ- 
ence, whereas positive peak values were found in the com- 
pressed domains. (A similar dependence of the layer align- 
ment in flow direction depending on the system size has been 
recently reported in^^.) Smaller systems with 16 layers across 
the gradient direction, which is comparable with om- systems 
consisting of A'^, x Ny x N; = 256 x 128 x 1 sites, showed a 
larger tendency to align in flow direction. For larger systems 
there was a clear trend to align at an angle with the flow di- 
rection. These distortions of the order parameter contribute a 
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term in f to the shear stress whose magnitude could be di- 
rectly related to the density of defects. In one system we found 
shear thinning described by a power law dependence of stress 
on strain rate, = Y" with m = 0.87; however the simplicity 
of this result may be fortuitous since the range of strain rates 
studied spans %. In another system, at strain rates primarily 
below %, the shear rheology was quasi-Newtonian (m = 1) al- 
though a significant first normal stress difference A^i was also 
present. 

We then simulated the flow behaviour of directly compa- 
rable systems in three dimensions, and found striking differ- 
ences between the two- and three-dimensional response under 
conditions of matched y. The 3D-system showed none of the 
distictive dilation/compression bands seen in 2D; moreover 
the microstructure, initially amorphous after quenching from 
the homogeneous phase, rapidly ordered under shear into a 
lamellar array with a low density of edge defects, all of which 
were oriented along the flow directions. This contrasts with 
the 2D case where numerous edge dislocations are present in 
the disordered phase (y > y^ ) but, by virtue of the 2D geom- 
etry, all are oriented along the vorticity direction. (In saying 
this, we imagine a 3D continuation of the 2D structure which 
is translationally invariant along that direction.) Our 3D re- 
sults do not rule out a transition or crossover from ordered to 
disordered structures, but this is not seen within the range of 
flow rates addressed here. This means that, if such a transi- 
tion is controlled by capillary number, its value at onset in 3D 
obeys Ca^ > 0.448, a value notably higher than that observed 
in 2D. Since identical parameter sets were chosen, the same 
applies to any other dimensionless number that might control 
the onset of the defect-rich phase. Finally, we note that one 
recent theory of smectic rheology suggests this is controlled 
by screw dislocations which have no 2D analogue. Moreover, 
the residual dislocations seen in the 3D ordered phase have an 
orientation that also does not exist in 2D. It should be remem- 
bered that, in 3D, edge dislocations and screw dislocations are 
interconvertible: they have the same topological status and the 
nomenclature merely describes whether the defect line is ori- 
ented normal or parallel to the smectic layers'^. The ability 
of 3D defect lines to reorient as well as translate underscores 
their greater dynamical freedom compared to similar defects 
in a 2D geometry. It is thus tempting to attribute the marked 
difference in dynamics between 2D and 3D systems to the re- 
duced constraints on defect motion that arise in 3D, allowing 
them to find and annihilate each other more easily. 
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7 Appendix: Sliding Periodic Boundary Con- 
ditions 

The simulations of shear flow are undertaken using a mixed fi- 
nite difference / lattice Boltzmann approach, where the former 
is used to advance the order parameter in time via the Cahn- 
Hilliard equation, and the latter is used to solve the Navier- 
Stokes equations. Coupling between the two is via the veloc- 
ity field, which advects the order parameter, and the thermody- 
namic pressure tensor (6), the divergence of which contributes 
a body force on the fluid. The sliding periodic boundary con- 
ditions of Lees and Edwards'*^- rnaybe applied to the lattice 
Boltzmann part of the calculationl^SSI to generate shear in the 
system. We describe here how the sliding periodic boundaries 
are extended to include the finite difference part of the calcu- 
lation. 

We consider a finite difference mesh where the order pa- 
rameter values are co-located with the density and velocity 
fields provided by lattice Boltzmann. Away from the sliding 
periodic boundaries, the finite difference discretisation is as 
one would expect. In order to ensure conservation of order 
parameter, as expressed by 

M + Mua<i?-Mdap)=0, (13) 

we consider each lattice point to be surrounded by a cubic 
control volume or cell of dimensions Ax = 1 in two or three 
dimensions. We compute convective and diffusive fluxes at 
the cell boundaries using linear interpolation of the normal 
velocity component to the cell face. For the advective fluxes, 
the interfacial values of the order parameter are computed via 
cubic interpolation weighted in the upwind direction. This is 
combined with an Euler forward step for the time integration: 
a small Courant number (uAt/Ax), and hence numerically sta- 
bility, is guaranteed by the small Mach number constraint stan- 
dard in lattice Boltzmann. The gradient of the chemical poten- 
tial is computed at each cell face using values of /j computed 
at the lattice point either side. 

A body force on the fluid in the Navier-Stokes equations 
may be included in the lattice Boltzmann approach. Again, to 
ensure conservation of momentum, the divergence of the ther- 
modynamic pressure tensor is computed based on cell-face 
values. In practice, this requires a finite difference stencil that 
entends at least ±3 lattice points in each coordinate direction 
(the thermodynamic pressure involves a fourth derivative of 
the order parameter). 

The sliding periodic boundaries lie (conceptually) along 
plane co-incident with the cell faces. Where the finite dif- 
ference stencil overlaps a such a plane, the relative movement 
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between the two sides of the plane must be accounted for when 
computing derivatives of the order parameter, and hence fluxes 
of order parameter and momentum, near the planes. This 
means interpolating values of the order parameter on each side 
of the plane to positions corresponding to the lattice locations 
on the other side so that the finite difference stencil can be 
used normally. For this study, where there are high deriva- 
tives of the order parameter in the free energy, it is important 
that this interpolation is at least cubic. Linear interpolation of 
the order parameter leads to unphysical oscillations of physi- 
cal quantities such as the shear stress related to the magnitude 
of the velocity jump across the sliding boundary. With cubic 
interpolation, such artefacts are absent. The non-linear nature 
of this interpolation means that fluxes of both order parameter 
and momentum computed at the plane do not exactly match 
on each side. However, this may be corrected at each sliding 
boundary at each time step, restoring global order parameter 
and momentum conservation. 
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